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Trapped, laser-cooled atoms and ions are quantum systems which can be experimentally con- 
trolled with an as yet unmatched degree of precision. Due to the control of the motion and the 
internal degrees of freedom, these quantum systems can be adequately described by a well known 
Hamiltonian. In this colloquium, we present powerful numerical tools for the optimization of the 
external control of the motional and internal states of trapped neutral atoms, explicitly applied to 
the case of trapped laser-cooled ions in a segmented ion-trap. We then delve into solving inverse 
problems, when optimizing trapping potentials for ions. Our presentation is complemented by a 
quantum mechanical treatment of the wavepacket dynamics of a trapped ion. Efficient numer- 
ical solvers for both time-independent and time-dependent problems are provided. Shaping the 
motional wavcfunctions and optimizing a quantum gate is realized by the application of quantum 
optimal control techniques. The numerical methods presented can also be used to gain an intu- 
itive understanding of quantum experiments with trapped ions by performing virtual simulated 
experiments on a personal computer. Code and executables are supplied as supplementary online 
materia^ 
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I. INTRODUCTION 

The information carrier used in computers is a bit, 
representing a binary state of either zero or one. In the 
quantum world a two-level system can be in any super- 
position of the ground and the excited state. The basic 
information carrier encoded by such a system is called a 
quantum bit (qubit). Qubits are manipulated by quan- 
tum gates — unitary transformations on single qubits and 
on pairs of qubits. 

Trapped ions are among the most promising physical 
systems for implementing quantum computation ( |Nielsen| 
and Chuang 2000 ) . Long coherence times and individual 



addressing allow for the experimental implementation of 
quantum gates and quantum c omputing protocols such 
as the Deu tsch-Josza algorithm ( Guide et al. 2003 1 , tele- 
portation ( Barrett et al. 2004 Riebe et al\ 2004 ) , quan- 



tum error correction (Chiaverini et al. 2004), quantum 



Fourier transform (Chiaverini et al. 20051 and Graver's 



search (Brickman et al. 2005p . Complementary research 



is using trapped neutral atoms in micro potentials such as 



2009 



2001 



magnetic micro traps (Fortagh and Zimmermann 2007 



Schmiedmayer et al. 20021, dipole traps flGaetan et al. 



Grimm et al. 2000p o r optical lattices ( Lukin et al. 



Mandel et al. 2003) where even individual imag- 



ing of single atoms has b een accomplished ( |Bakr et al. 



2009 Nelson et al. 2007). The current challenge for all 



approaches is to scale the technology up for a larger num- 



2 



ber of qubits, for which several proposals exist ( 


Cirac and 


Zoller 


2000 Duan et al. 


2004 


Kielpinski et al. 


2004 


)• 



The basic principle for quantum computation with 
trapped ions is to use the internal electronic states of 
the ion as the qubit carrier. Computational operations 
can then be performed by manipulating the ions by co- 



herent laser light (Blatt and Wineland 2008 Haffner 



et al. 2008). In order to perform entangling gates be- 



tween different ions, Cirac and Zoller ( 1995 ) proposed 



to use the mutual Coulomb interaction to realize a col- 
lective quantum bus (a term used to denote an object 
that can transfer quantum information between subsys- 
tems). The coupling between laser light and ion motion 
enables the coherent mapping of quantum information 
between internal and motional degrees of freedom of an 
ion chain. Two-ion gates are of particular importance 
since combined with single-qubit rotations, they consti- 
tute a universal se t of quantum gates for computation 
(DiVincenzo 1995). Several gate realizations have been 



prop osed (|Cirac and Zoller 



2005 



Jonathan et al. 



and Wunderlich, 2001 



2000 tv 



1995 



Garci'a-Ripoll et al. 



ilburn et al. 1 12000; Mintert 



M0fmer and S0rensen| [l999 Mon- 



roe et al. 1997 |Poyatos et al. 1998[) and realized b y sev 



eral groups (|Benhelm et al 
Kim et aT\ 



2008; DeMarco et al. 2002 



2008[ |Leibfried et aLfpOOBT 



Schmidt-Kaler 



et al. 2003c I . When more ions are added to the ion 
chain, the same procedure can be applied until the dif- 
ferent vibrational-mode frequencies become too close to 
be individually addressable 1 ; the current state-of-the-art 
is the preparation and read-out of an W entangled state 



of eight ions (Haffner et al. 2005) and a six- ion GHZ 



state (Leibfried et al. 2005) 



A way to solve this scalability problem is to use seg- 
mented ion traps consisting of different regions between 
which the ion s are shuttled (transported ba ck and forth) 
( |Amini et al] |2010| |Kielpinski et al.\ |2004[ ) . The proper 
optimization of the shuttling processes and optimization 
of the laser ion interaction can only be fully performed 



with the aid of numerical tools (Huber et al. 2010 Hu- 


cul et al. 


2008 


Reichle et al. 


2006[ Schulz et al. 


2006). 



In our presentation equal emphasis is put on the pre- 
sentation of the physics of quantum information exper- 
iments with ions and the basic ideas of the numerical 
methods. All tools are demonstrated with ion trap ex- 
periments, such that the reader can easily extend and 
apply the methods to other fields of physics. Included 
is supplementary material, e.g. source code and data 
such that even an inexperienced reader may apply the 
numerical tools and adjust them for his needs. While 
some readers might aim at understanding and learning 
the numerical methods by looking at our specific ion 
trap example others might intend to get a deeper under- 
standing of the physics of quantum information exper- 



1 although there are schemes that address multiple modes |Kim| 
et al\ 2009||Zhu et aL||2006a|b l 



iments through simulations and simulated experiments. 
We start in Sec. [TT] with the description of the ion trap 
principles and introduce numerical methods to solve for 
the electrostatic potentials arising from the trapping elec- 
trodes. Accurate potentials are needed to numerically 
integrate the equation of motion of ions inside the trap. 
Efficient stable solvers are presented in Sec. |III| The axial 
motion of the ion is controlled by changing the dc volt- 
ages of the electrodes. However, usually we would like 
to perform the inverse, such that we find the voltages 
needed to be applied to the electrodes in order to produce 
a certain shape of the potential to place the ion at a spe- 
cific position with the desired trap frequency as described 
in Sec. |IV| This problem belongs to a type of problems 
known as inverse problems, which are quite common in 
physics. In Sec. [V] we enter the quantum world where 
we first will obtain the stationary motional eigenstates 
of the time-independent Schrodinger equation in arbi- 
trary potentials. We then describe methods to tackle 
the time-dependent problem, and present efficient nu- 
merical methods to solve the time-dependent Schrodinger 
equation. The presented methods are used in Sec. |VI| 
where we consider time-dependent electrostatic poten- 
tials with the goal to perform quantum control on the 
motional wavefunction and present the optimal control 
algorithm. Finally, we apply these techniques in Sec. |VII| 
to the Cirac-Zoller gate. In the conclusion Sec. |VIII[ we 
give a short account on the applicability of the presented 
numerical methods to qubit implementations other than 
trapped laser cooled ions. 



II. ION TRAP DEVELOPMENT 
ELECTROSTATIC FIELDS 



CALCULATION OF 



The scalability problem for quantum information with 
ion traps can be resolved with segmented ion traps. The 
trap potentials have to be tailored to control the posi- 
tion and trapping frequency of the ions. In the follow- 
ing, we will describe the mode of operation of a simple 
ion trap and then present numerical solvers for the effi- 
cient calculation of accurate electrostatic fields. Due to 
the impossibility of generating an electrostatic potential 
minimum in free space, ions may either be trapped in 
a combination of electric and magnetic sta tic fields - a 
Penning trap (Brown and Gabrielse 1986), or in a ra- 



dio frequency electric field - a Paul trap, where a radio 
frequency (rf) voltage U T [ with rf drive frequency tu r { is 



applied to some of the ion-trap electrodes (Paul 
In the latter case, we generate a potential 

$(x,y,z,t) = Lr ^{a dc x 2 + f3 dc y 2 +j dc z 2 ) 



1990). 



2 
2 



3(w r fi)(a rf x 2 + /3 rl y 2 + 7 rf z 2 ),(l) 

where U dc is a constant trapping voltages applied to 
the electrodes. The Laplace equation in free space 
A$(x,y, z) = puts an additional constraint on the co- 
efficients: a dc + p dc + 7dc = and a r f + /3 T t + j x f — 0. One 



possibility to fulfill these conditions is to set aye = Pdc = 
7d c = and a T t + /3 r f = — 7 r f. This produces a purely 
dynamic confinement of the ion and is realized by an 
electrode configuration as shown in Fig. [lja) , where the 
torus-shaped electrode is supplied with radio frequency 
and the spherical electrodes are grounded. An alterna- 
tive solution would be the choice — aye = /?dc + 7dc an d 
arf = 0, /3,f = — 7rf, leading to a linear Paul-trap with dc 
confinement along the x-axis and dynamic confinement 
in the yz-plane. Fig. [TJb) shows a possible setup with 
cylindrically shaped electrodes and segmented dc elec- 
trodes along the axial direction which we will consider 
in the following. In this trapping geometry, the ions can 
crystallize into linear ion strings aligned along the re-axis. 
The classical equation of motion for an ion with mass m 
and charge q is mx = — qV$, with x = (x, y, z) (James 
1998). For a potential given by Eq. ([IJ the classical equa- 



tions of motion are transformed into a set of two uncou- 



pled Mathieu differen tial equations ( |Haffner et al. 



Leibfried et al. 2003 1 



2008 



+ K - 2q u cos(20)ti(f) = u = y,z, (2) 



with 2£ = w r fi. The Mathieu equation belongs to the 
family of differential equations with periodic boundary 
conditions and its solution is readily found in textbooks 
(for example, (Abramowitz and Stegun 1964 1). For a 



linear Paul-trap, the parameters a v 
plane are given by 



and q u in the yz- 



Qz 



2\q\Urif3 lf 

2 ' 

2|g|£7rf7rf 



4|g|t/dc/?dc 
4Mt/dc7dc 



(3) 



The sol ution is sta ble in the range < /3 U < 1, where 
/?„ = \J a u + q^/2 only depends on the parameters a u 
and q u . The solution of Eq. ([2| in the lowest order ap- 
proximation (| an |, 9^ <C 1), which implies that j3 u <C 1, 
is 



u(t) — u cos(cj u £) (l + ^ cos(u) x ft)J 



(4) 



The ion undergoes harmonic oscillations at the secu- 
lar frequency uj u — /? u o» r f/2 modulated by small os- 
cillations near the rf-drive frequency (called micromo- 
tion). The static axial confinement along the a;-axis is 
harmonic with the oscillator frequency being given by 
uj x = ■\/|q|C/dcQ ! dc/ m - The axial confinement is gener- 
ated by biasing the dc electrode segments appropriately. 
Typical potential shapes can be seen in Fig. [2]ja). 

The radial confinement is dominated by the rf potential 
which can be approximated by an effective harmonic po- 
tential $ e &(y> z ) = M |V$(y,z)|V(4™w r 2 f ) where $(y,z) 
is the potential generated by setting the radio frequency 
electrodes to a constant voltage U T f see Fig. [2j». How- 
ever, this effective potential is only an approximation and 
does not take the full dynamics of the ion into account. 




FIG. 1 (Color online). Electrode geometries of ion traps: The 
rf electrodes are depicted in blue and dc electrodes in yellow 
respectively, (a) Typical electrode configuration for a 3D ring 
trap with dynamic rf confinement in all three dimensions, (b) 
Electrode arrangement for a linear Paul trap. The dc elec- 
trodes are divided into segments numbered from 1 to 5. For 
the numerical simulations we assume the following parame- 
ters: Segment have a width of 2 mm and a radius of 0.5 mm. 
The central dc electrode is centered at the position x = 0. 
The minimum distance of the electrode surface to the trap 
axis is 1.5 mm. 
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FIG. 2 (a) Trapping potentials along the a;-axis generated 
by each individual electrode from the linear Paul trap geom- 
etry of Fig. [TJb). Each curve corresponds to the respective 
electrode biased to -1 V and all others to V. (b) Equipoten- 
tial lines of the pseudo-potential in the radial plane (U r f =200 
V pp , w r f = 2-7T x 20 MHz). Potentials are obtained as described 
in Sec. HI] 



Before we can simulate the motion of the ion we need 
fast and accurate electrostatic field solvers. In the next 
section we first present the finite difference method and 
then the finite element method. If the potentials need 
to be known on a very small scale, a huge spatial grid 
would be needed. Therefore, we introduce the bound- 
ary element method and show how the efficiency of this 
method can be drastically improved by the application 
of the fast multipole method. 



A. Finite difference method 

To obtain the electrostatic potential &(x,y,z) in free 
space generated by a specific voltage configuration Ui 
for i = 1, ...,n applied to the n electrodes, we need 
to solve the Laplace equation A.G>(x, y, z) = 0, with the 



4 



Dirichlet boundary condition <&(x,y,z) = f7, for points 
lying on the ith electrode. There are several approaches 
to obtain the solution. The most intuitive is the finite 
difference method (FDM). The principle is that we can 
write the differential equation in terms of finite differ- 
To illustrate this, we take the one 



ences (Thomas 



1995) 



dimensional differential equation ^ = F(x) with the 
boundary condition $(0) = a where F (x) is an arbitrary 
function. If we write 



dx 



lim 

Ax-tO 



$(x + Ax) - $(x) 



F(x), 



(5) 



take only a finite difference Ax and discrctizc the x-axis 
by defining Xi — i Ax with i running from to N and 
xn = 1, we obtain a discrete approximation which di- 
rectly gives an explicit update equation (using the Euler 
method) for $: 



$(x l+ i) = $(xi) + AxF(x,)- 



(6) 



Eq. ^ can then be applied iteratively to solve the differ- 
ential equation. By comparing the solution with the Tay- 
lor expansion and assuming that the higher order terms 
are bounded, we see that the error of this finite difference 
approximation is of order Ax. This is usually written as 



$(x + Ax) - $(x) _ d$ 
Ax dx 



O(Ax), 



(7) 



which means that there exists a constant d such that 



Ax 



d'I> 
dx 



< d|Ax| for all x. The Laplace 
equation is of second order, but one can transform it into 
a set of a first order differential equations 



equations. We have to specify two boundary conditions 
which we assume to be 3>(xn) = U\ and $(xjv) = Ui, 
where U\ and U2 are the voltages supplied at the bound- 
aries. The matrix equation then has the form 



/-2 1 ■■• 0\ 

1-2 1 : 
1 -2 '•■ 



■• 1 

1 -2/ 



/ $ (^i) \ 
*(as 2 ) 
*(as 3 ) 

Wx N ^)J 





\-uJ 



(12) 



This equation has a tridiagonal form and can be most 



efficie ntly solved by the Thomas algorithm (Press et al. 
2007 ) 2 . A sparse matrix with more off-diagonal entries 



is obtained when the two- or three-dimensional Laplace 
equation is treated in a similar fashion. The solution is 
then obtained either by simple Gaussian elimination or 
more efficiently by iterative methods, such as the suc- 



cessive over relaxation method (SOR) ( Press et al. 2007 ) 
or the generalized minimum residual method (GMRES) 
dSaadl|2bol. 



One advantage of FDM is that it is easy to implement 
on a uniform Cartesian grid. But in modeling three- 
dimensional geometries one usually favors a triangular 
non-uniform mesh, where the mesh spacing is spatially 
adapted to the local complexity of the geometry struc- 
tures; e.g. it makes sense to use a finer mesh near the 
edges. 



d /$ 
dx \v 



v 

F(x) 



(8) 



from which an explicit update rule can be derived, which 
is O(Ax). We can obtain a second order approximation 
by cancelling the first order terms which gives a centered- 
difference approximation for the first derivative 



$(x 



n+l) 



$(X„_! 



2Ax 



dx 



C(Ax 2 



(9) 



which is of order C(Ax 2 ). A centered-difference approx- 
imation for the second derivative reads 



$(x» + i) - 2$(x») + $(x»_i) 
Ax 2 



d 2 $ 
dx 2 



0(Ax 2 ), 



(10) 

which is again of order C(Ax 2 ). The update rule for the 
one-dimensional Laplace equation 4-^ = thus has the 



form 



$(x„ +1 ) - 2$(x„) + $(x n _i) = 0, 



(11) 



B. Finite element method 

The finite element method (FEM) is better suited 
for nonuniform meshes with inhomogeneous granularity, 
since it transforms the differential equation into an equiv- 
alent variational one: instead of approximating the differ- 
ential equation by a finite difference, the FEM solution 
is approximated by a finite linear combination of basis 
functions. Again, we demonstrate the method with a 
one-dimensional differential equation = F(x), and 
for simplicity we take the boundary condition $(0) = 
and $(1) = 0. The variational equivalent is an integral 
equation integrated between the boundaries at and 1: 



d 2 $(x) 
dx 2 



v{x)dx 



F{x)v(x)dx, 



(13) 



where v(x) is the variational function which can be freely 
chosen except for the requirement v(0) = v(l) = 0. Inte- 



which is an implicit expression, since the solution now 
has to be obtained by solving a linear system of algebraic 



2 see package octtool, function tridag 



5 




FIG. 3 (Color online). Overlapping basis functions from 
Eq. ( |15[ ) Vk(x) (colored solid lines) for the finite element 
method providing linear interpolation (black dashed line) of 
an arbitrary function (black solid line). 



grating this by parts gives 



i A2 



d 2 >S>{x) 
dx 2 



v(x)dx 



dx 



v(x) 



1 d$(x) dv(x) 



dx dx 



dx 



F(x)v(x)dx. 



(14) 



We can now discretize this equation by constructing v(x) 
on a finite-dimensional basis. One possibility is linear 
interpolation: 



v k {x) 



X k —X k -1 



Xk-1 < X < Xk, 
X k < X < X k+ 1, 



otherwise 



(15) 



with xq — 0, xn = 1 and Xk are the (not neces- 
sarily equidistant) sequential points in between and k 
ranges from 1 to N — 1. These functions are shown 
in Fig. [3j The advantage of this choice is that the in- 
ner products of the basis functions f Vk(x)vj(x)dx and 

their derivatives f Q v' k {x)v'Ax)dx are only nonzero for 
\j — k\ < 1. The function and F(x) are then ap- 

proximated by <I>(:r) w SfcLo &( x k)vk(x) and F(x) w 
SfeLo F{ x k)vk(x) which linearly interpolates the initial 
functions (see Fig. §. With ^ » £f 
Eq. ( 14 ) is now recast into the form 



N 



k=0 



dvk (x) dvj (x) 
dx dx 



N r f 1 

F(x k ) / v k {x)v J (x)dx 
k=o L Jo 



(16) 



where the terms in brackets are sparse matrices. This 
matrix equation can then again be solved by iterative 
matrix solvers (such as GMRES). 

For the Laplace problem, we need to extend this 
method to higher dimensions. In this case, instead of the 



integration by parts in Eq. (14) we have to use Green's 
theorem (Jackson 20091: 



<9$ 

A$(x)v(x)dV = I —vds- 

Igy on 



F(x)v(x)dV, 



(17) 



where V is the volume of interest and SV the bound- 
ing surface of the volume. Now space is discretized by 
three dimensional basis functions and we can proceed in 
an analogous manner as in the one dimensional case de- 
scribed above. 

Potentials obtained by FDM and FEM usually result 
in unphysical discontinuities (i.e. numerical artifacts) and 
must be smoothed in order to be useful for ion trajectory 
simulations. Additionally, in order to obtain high accu- 
racy trajectory simulations needed to simulate the trajec- 
tory extend of a trapped ion of less than 100 nm, the po- 
tentials that are calculated have to be interpolated, since 
computing with a grid with nanometer spacing would in- 
volve an unbearable computational overhead: the whole 
space including the typically centimeter sized trap would 
have to be meshed with a nanometer-spaced grid. FEM 
would allow for a finer mesh in the region where the ion 
would be located reducing the overhead somewhat, but 
this does not increase the accuracy of the surrounding 
coarser grid. Avoiding to give a wrong expression we 
would like to stress that the FEM method finds wide 
applications in engineering and physics especially when 
complicated boundary conditions are imposed but for our 
accuracy goals FEM and FDM are inadequate. 



C. Boundary element method - fast multipole method 

We proceed to show a different way of solving the 
Laplace problem with a method which features a high 
accuracy and gives smooth potentials that perform well 
in high-resolution ion-ray-tracing simulations. 

To begin with, we divide the electrodes into small sur- 
face elements Si of uniform surface charge density <7j, 
with i numbering all surface elements from 1 to N. The 
potential at any point in space caused by a charge dis- 
tribution of these elements can be easily obtained from 
Coulomb's law: one must simply sum up all the contri- 
butions from each surface element. Hence the voltage Uj 
on the surface element Sj is generated by a linear super- 
position of the surface charge densities (Ji — d$(xi)/dn 
also expressed as the normal derivative of the potential <& 
(as obtained from the Maxwell equations) on all surface 
elements (including Sj) additionally weighted by geome- 
try factors given by the Coulomb law (represented by a 



G 



matrix C) providing the following simple matrix equation 



(18) 



Now we want to solve for the surface charge densities in 
terms of the applied voltages. The surface charge densi- 
ties for a given voltage configuration can then be obtained 
by finding the matrix inversion of C . This is the basic 



Fortunately, Greengard and Rokhlin ( 1988 ) came up 

















= c 




\u N ) 




{ctnJ 



idea of the boundary element method (BEM) ( Pozrikidis 



2002). In the case of metallic surface elements where ei- 
ther the potential or the charge density is fixed, we have 
to exploit Green's second identity 



JY 



rate 

On 



N 

-2^ft(x J )<f(x i ). (19) 

i=l 



This equation has then to be solved for the unknown pa- 
rameters which are the surface charge density 3 ^ c '- ) on 
surfaces with given potential or the potential $(xj) on 
surfaces with given charge density. Now we can choose 
the Xj and Xj to be representative points on each surface 
element e.g. the center of gravity. This corresponds to 
the approximation that the potential and charge density 
are constant on each surface element. Eq. (19) is a ma- 
trix equation equivalent to Eq. (18). Oii is obtained by 



performing a surface integral over the surface elements 
Si of the two-dimensional Green's function G(x,Xj) = 
— ^- In |x — Xj| (for two-dimensional problems) or the 
three-dimensional Green's function G(x, Xj) = 47r |~^ x .| 
(for three-dimensional problems), ft is obtained by per- 
forming a surface integral over the surface elements Sj 
of the gradient of the Green's function multiplied by the 
surface norm n 



*i(xj) = <j> G(x,Xj)ds, 
<j> n(x) ■ VG(x,Xj-)ds. 



(20) 



(21) 



Analytical expressions for these integrals over triangu- 



lar surface elements can be found in |Davey and Hinduja| 
( 1989 ) , or via Gauss-Legendre quadrature over a trian- 



gle. Eq. ( 19 ) is now solved for the unknown parameters 
such as the surface charge densities . Once this is 

achieved, we can calculate the potential 



$(x) = -^a j (x)-A-^ 



JV 



+ > Mx)*(xi) (22) 



t=i 



at any position x with ccj(x) and ft(x) evaluated at the 
same position. BEM is very accurate and the implemen- 
tation is quite straight-forward, but the complexity of the 
matrix inversion scales prohibitively as 0(N 3 ). Different 
to the finite element method we cannot use sparse matrix 
solvers for this matrix inversion. 



with an innovative method for speeding-up the ma- 
trix vector multiplication needed for iterative matrix in- 
version, which they termed the fast multipole method 
(FMM). FMM can solve the BEM problem with 0(N) 
complexity, giving a drastic increase in speed, and mak- 
ing BEM applicable to more complex systems. In a se- 
ries of publications, the algorithm was further improved 



( Carrier et al. 


1988 


1999 Greengard and Rokhlin 


1997 


Gumcrov anc 


Duraiswami[ |2005| |Nabors et al.\ 


1994 


Shen and Liu 


2007 


) and extended to work with the 


Helmholtz equation 


(Gumerov and Duraiswami 2004|). 



The basic idea was to use local and far field multipole 
expansions together with efficient translation operations 
to calculate approximations of the fields where the three- 
dimensional space is recursively subdivided into cubes. A 
detailed description of the method is beyond the scope 
of this paper and we refer to the cited literature. 



D. Application 



We have used the FMM implementation from |Nabors| 
et al. ( 1994[ ) and combined it with a scripting language for 
geometry description and the ability to read AutoCAD 
files for importing geometrical structures. Any small in- 
accuracies due to numerical noise on the surface charges 
are 'blurred out' at large distances due to the Coulomb 
law's 1/r scaling. In this regard, we can assert that the 
surface charge densities obtained by FMM are accurate 
enough for our purposes. If special symmetry proper- 
ties are needed (such as rotational symmetry for ion-lens 
systems or mirror symmetry) then one can additionally 
symmetrize the surface charge densities. We have imple- 
mented symmetriz ation functions in o ur code to support 

As FMM is used 



these calculations (Fickler et al 



2009) 



to speed up the matrix vector multiplication it can be also 
used to speed up the evaluation of Eq. ( p2| to obtain the 
potentials in free space. However, if accurate potentials 
in the sub micrometer scale are needed (such as for our 
application), it is better to use FMM for the calculation 
of the surface charge densities i. e. for the inversion of ma- 
trix Eq. ( 19 ) and then use conventional matrix multipli- 
cation for the field evaluations as described by Eq. (22). 
Fig.[2]ja) shows the smooth potentials calculated by solv- 
ing for the surface charge densities with FMM. Depicted 
are the potentials for each electrode when biased to -1 V 
with all others grounded. A trapping potential is then 
generated by taking a linear superposition of these po- 
tentials. Fig. |2jb) shows the equipotential lines of the 
pseudo-potential. The full implementation can be found 
inside our bemsolver package together with example files 
for different trap geometries. 

With the calculated potentials from this chapter we 
can now solve for the motion of an ion in the dynamic 
trapping potential of the Paul trap which will be the 
focus of the next chapter. 
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III. ION TRAJECTORIES 
MOTION 



CLASSICAL EQUATIONS OF 



The electrostatic potentials obtained with the methods 
presented in the previous chapter are used in this section 
to simulate the trajectories of ions inside a dynamic trap- 
ping potential of a linear Paul trap. We present the Eulcr 
method and the more accurate Runge-Kutta integrators. 
Then we show that the accuracy of trajectories can be 
greatly enhanced by using phase space area conserving 
and energy conserving solvers such as the Stormer- Verlet 
method, which is a partitioned Runge-Kutta integrator. 



A. Euler method 

The equation of motion of a charged particle with 
charge q and mass m in an external electrical field can 
be obtained by solving the ordinary differential equation 



x(t) = f(t,x) 

x\ / V 



y(*) = 



f(x,i) 



where f(i,x) = (q/m)E(t, x) = -(g/m)V$(t,x) is the 
force arising from the electric field. The vectors y and 
F are six-dimensional vectors containing the phase space 
coordinates. As in the previous section, the equation 
of motion can be solved by means of the explicit Eulcr 
method with the update rule y n +i = Yn + hF(t n ,y n ), 
where we use the notation y„ = y(i„). If s is the ab- 
solute tolerable error, then the time step h = t n+1 — t n 
should be chosen as h = yfs, which gives the best com- 
promise between numerical errors caused by the method 
and floating point errors accumulated by all iterations. 

An implicit variation of the Euler method is given 
by the update rule y n+ i = y„ + h¥(t n+u y n+1 ). Nei- 
ther of the methods is symmetric, which means that 
under time inversion (h — » —h and y„ — > y n +i), a 
slightly different trajectory is generated. A symmet- 
ric update rule is given by the implicit midpoint rule 
Yn+i = y« + hF((t n + i„+i)/2, (y„ + y„ + i)/2), which 
has the additional property that it is symplectic, mean- 
ing that it is area preserving in phase space. The explicit 
and implicit Euler methods are of order O(h) whereas 
the implicit midpoint rule is of order 0(h 2 ) (Hairer et al. 
2002). These methods belong to the cla ss of one stage 
Runge-Kutta methods ( Greenspan] 2006). 



B. Runge-Kutta method 

The general s-stage Runge-Kutta method is defined by 
the update equation 



y«+i = y n + h ^2 & » k * 

1=1 



(23) 
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TABLE I Butcher tableau for the 4th and 5th order Runge- 
Kutta methods: The left most column contains the d coeffi- 
cients, the last two rows under the separation line contain the 
h coefficients to realize a 4th order or 5th order Dormand- 
Price Runge-Kutta. The atj coefficients are given by the re- 
maining numbers in the central region of the tableau. Empty 
entries correspond to = 0. 



with 



k, = F(t n 



h (Hj 

3=1 



3=1 



■'J ' 



(24) 



for i — 1 , . . . , s and hi and real numbers, which 

are given in Butcher tableaux for several Runge-Kutta 
methods. Note that in the general case the ki are defined 
implicitly such that Eq. (24) have to be solved at each 
time step. However, if a,j = for i < j then the Runge- 
Kutta method is known as explicit. The standard solver 
used in many numerical packages is the explicit 4th and 
5th order Dormand- Price Runge-Kutta, whose values are 
given in the Butcher tableau in Tab. [TJ The difference 
between the 4th and the 5th order terms can be used as 
error estimate for dynamic step size adjustment. 



C. Partitioned Runge-Kutta method 

A significant improvement can be achieved by parti- 
tioning the dynamical variables into two groups e.g. po- 
sition and velocity coordinates, and to use two different 
Runge-Kutta methods for their propagation. To illus- 
trate this we are dealing with a classical non-relativistic 
particle of mass m which can be described by the fol- 
lowing Hamilton function iJ(x, v) = T(v) + $(x) with 
T(v) = mv 2 /2. The finite difference version of the 
equation of motion generated by this Hamiltonian reads 
x n _|_i — 2x„ + X„_ i = h 2 { (x„). The only problem is that 
we cannot start this iteration, as we do not know x_i. 
The solution is to introduce the velocity v = x written 
as a symmetric finite difference 



2h 



(25) 



and the initial conditions: x(0) = x , x(0) = v such 
that we can eliminate x_i to obtain xi = xq + /ivq + 
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TABLE II Butcher tableau for the 4th order partitioned 
Runge-Kutta method consisting of a 3 stage Lobatto IIIA- 
IIIB pair (Hai rer et a/.| [2 002 |: Left table shows the 6j, 
and Ci coefficients corresponding to the table entries as de- 
scribed in the legend of Tab.|l Right table shows correspond- 
ing primed b[, a[j and c' t coefficients. 



4-f (x ). Now we can use the following recursion relation 

h 



v„+i/2 = v„ + -f(x„), 



Wi+l 



hv. 



+ 1/2! 



V n+ 1 = V n+ i/ 2 + -f(x„ + i). 



(26) 
(27) 

(28) 



One should not be confused by the occurrence of half in- 



teger intermediate time steps. Eqs. (26) and (28) can 
be merged to v n+1 / 2 = v„_ 1/2 + /if(xT) if v n+ i is not 
of interest. The described method is the Stormer-Verlet 
method and is very popular in molecular dynamics sim- 
ulations where the Hamiltonian has the required proper- 
ties, e.g. for a conservative system of iV particles with 
two-body interactions: 



ff(x,v) 



^ N N i—1 



2 ^ 



(29) 



i=2 j=l 



In the case of the simulation of ion crystals, Vij would be 
the Coulomb interaction between ion pairs. The popular- 
ity of this method is due to the fact that it respects the 
conservation of energy and is a symmetric and symplectic 
solver of order two. It belongs to the general class of par- 
titioned Runge-Kutta methods where the s-stage method 
for the partitioned differential equation y(i) = f(t, y, v) 
and v(t) = g(t, y, v) is defined by 



Yn+i = y„ + hy^b t kj, 

i=i 

s 

v„+i = v n + hy^b'jlj, 



(30) 
(31) 



kj = f (t n + Oj h, y n + h Ojj kj , v n + h } j o< 7 - lj ) , 



3=1 



3 = 1 



3=1 3 = 1 

with bi and atj = 1, . . . , s) being real numbers and 
= X)j=i a ij an d analogous definitions for b^, a'^ and c^. 
As an example Tab. ITT] shows the Butcher tableau for the 
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FIG. 4 Comparison of the Euler and Stormer-Verlet simula- 
tion methods: (a) Shows trajectories in the radial plane sim- 
ulated with the explicit Euler method. It can be clearly seen 
that the trajectories are unstable, (b) Simulation for the same 
parameters with the Stormer-Verlet method results in stable 
trajectories. The small oscillations are due to the micro mo- 
tion caused by the rf drive, (c) Phase space trajectories of the 
harmonic axial motion of an ion numerically integrated with 
the explicit Euler method. The simulation was performed for 
a set of three different initial phase space coordinates, as indi- 
cated by the triangles. Energy and phase space area conser- 
vation are violated, (d) Equivalent trajectories numerically 
integrated with the Stormer-Verlet method which respects 
energy and phase space area conservation. The parameters 
are: Simulation time 80 pis, number of simulation steps 4000, 
Urf =400 V pp , w rf = 2tt x 12 MHz, U dc = (0, 1,0, 1,0)V for 
the trap geometry shown in Fig. [JJb). 



4th order partitioned Runge-Kutta method consisting of 
a 3 stage Lobatto IIIA-IIIB pair. A collection of more 



Butcher tableaux can be found in (Hairer et al. 2002) 



and (Greenspan 2006) 



D. Application 

We have simulated the phase space trajectories of ions 
in the linear Paul trap of Fig.[TJb) for 500,000 steps with 
the Euler method and with the Stormer-Verlet method. 
A trajectory in the radial plane obtained with the Eu- 
ler method can be seen in Fig. |4^a), it clearly shows an 
unphysical instability. When simulated by the Stormer- 
Verlet method, the trajectories are stable, as can be seen 
in Fig.^b). The small oscillations are due to the micro 
motion caused by the rf drive. The Euler method does 
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FIG. 5 (Color online). Illustration of the regularization tech- 
nique: Suppression of the divergence at zero (red curve) of the 
1/s term. The black curve shows the Tikhonov regularization 
term, the singular behavior of the 1/s inverse is avoided and 
diverging values are replaced by values near zero. The blue 
curve tends towards one for s —> where 1/s is diverging. It 
is used to force diverging inverses towards a given value (see 
text). 



not lead to results obeying to energy and phase space 
area conservation (see Fig. [4jc) ) , whereas the Stormer- 
Verlet method conserves both quantities (see Fig. ptd)). 
The Euler method should be avoided when more than a 
few simulation steps are performed. The Stormer-Verlet 
integrator is implemented in the bemsolver package. 

In the next section we will find out how we can control 
the position and motion of the ion. 



space, which could be chosen, for example, on the trap 
axis. We would like to position the ion at a specific loca- 
tion in a harmonic potential with a desired curvature, i.e. 
trap frequency. This is a matrix inversion problem, since 
we have specified $j over the region of interest but we 
need to find the voltages Uj. The problem here is that 
M is much larger than N, such that the matrix Aij is 
over determined, and (due to some unrealizable features 
or numerical artifacts) the desired potential might 
not lie in the solution space. Hence a usual matrix in- 
version, in most cases, will give divergent results due to 
singularities in the inverse of . 

This class of problems is called inverse problems, and 
if singularities occur then they are called ill-conditioned 
inverse problems. 

In our case, we wish to determine the electrode voltages 
for a potential-well moving along the trap axis, there- 
fore a series of inverse problems have to be solved. As 
an additional constraint we require that the electrode- 
voltage update at each transport step is limited. In the 
following, we will describe how the Tikhonov regulariza- 
tion method can be employed for finding a matrix inver- 
sion av oiding singularities and fulfilling additional con - 
straints( Press et al. 2007 Tikhonov and Arsenin, 19771. 



A. Thikonov regularization 

For notational simplicity, we will always assume that 
A is a M x N dimensional matrix, the potentials along 
the trap axis are given by the vector 



$ = ($(xi),$(X2) 



(33) 



IV. TRANSPORT OPERATIONS WITH IONS - 
ILL-CONDITIONED INVERSE PROBLEMS 

If we wish to transport ions in multi-electrode geome- 
tries, the electrostatic potential has to be controlled. An 
important experimental constraint is that the applicable 
voltage range is always limited. Additionally, for the dy- 
namic change of the potentials, voltage sources are gen- 
erally band limited and therefore we need the voltages of 
the electrodes to change smoothly throughout the trans- 
port process. The starting point for the solution of this 
problem is to consider A(xi,j) — A^, a unitless matrix 
representing the potential on all points Xi when each elec- 
trode j is biased to IV whereas all other electrodes are 
kept at 0V (see Fig. |2ja)). Hence we can calculate the 
generated total potential $(2^) = at any position Xi 
by the linear superposition 

N 

•I', X- 1 '.' - i = h---,M, (32) 
i=i 

with N denoting the number of separately controllable 
electrodes and M being the number of grid points in 



and u = (U~i,U2, ■ ■ ■ ,Un) t is a vector containing the 
electrode voltages. Instead of solving the matrix equa- 
tion Au = the Thikonov method minimizes the resid- 
ual ||Au — <J?|| 2 with respect to the Euclidean norm. This 
alone could still lead to diverging values for some com- 
ponents of u which is cured by imposing an additional 
minimization constraint through the addition of a regu- 
larization term such as 

a||u|| 2 , (34) 

which penalizes diverging values. The larger the real 
valued weighting parameter a is chosen, the more this 
penalty is weighted and large values in u are suppressed 
at the expense that the residual might increase. Instead 
of resorting to numerical iterative minimizers a faster and 
more deterministic method is to perform a singular-value 
decomposition which decomposes the M x N matrix A 
into a product of three matrices A = USV T , where U and 
V are unitary matrices of dimension M x M and N x N, 
and S is a diagonal (however not quadratic) matrix with 
diagonal entries Si and dimension M x N. Singular-value 
decomposition routines are contained in many numerical 
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libraries, for example lapack 3 . The inverse is then given 
by the A' 1 = VS'U T where S' = S' 1 . The diagonal en- 
tries of 5" are related to those of S by s • = 1/ Si such that 
the singular terms can now be directly identified. The ad- 
vantage of the singular value decomposition now becomes 
clear: all the effects of singularities are contained only in 
the diagonal matrix S' . By addressing these singulari- 
ties (i.e., when Si — > 0) in the proper way, we avoid any 
divergent behavior in the voltages Uj. The easiest way 
to deal with the singularities would be to set all terms 
s'i above a certain threshold to zero. Thikonov however 
uses the smooth truncation function 



(35) 



(see Fig. [5] black curve) which behaves like 1/sj for large 
Si but tends to zero for vanishing Sj , providing a gradual 
cutoff. The truncation parameter a has the same mean- 
ing as above in the regularization term: the larger a the 
more the diverging values are forced towards zero, and 
if a = then the exact inverse will be calculated and 
diverging values are not suppressed at all. The required 
voltages are now obtained by 



u = VS'U T $. 



(36) 



These voltages fulfill the requirement to lie within some 
given technologically accessible voltage range, which can 
be attained by iteratively adjusting a. 

In the remainder of this section, we present an exten- 
sion of this method which is better adapted to our specific 
task of smoothly shuttling an ion between two trapping 
sites: instead of generally minimizing the electrode volt- 
ages, we would rather like to limit the changes in the volt- 
age with respect to the voltage configuration Uo which is 
already applied prior to a single shuttling step. Therefore 
the penalty function of Eq. (34 1 is to be replaced by 



a u- u 



(37) 



The application of this penalty is achieved through an 
additional term in Eq. ( 36 ) 
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FIG. 6 (Color online). Voltage configurations to put the 
minimum of a harmonic potential with fixed curvature into 
a different positions. The insets show the resulting potentials 
obtained by linear superposition of the individual electrode 
potentials. 



to transform the vector Uo into the basis of the matrix 
D, 

A remaining problem now arises when e.g. a single 
electrode, i.e. a column in Eq. (32) leads to a singular- 
ity. This might due to the fact that it does generate only 
small or vanishing fields at a trap site of interest. Even 



though the singularity is suppressed by the regularization 
it nevertheless leads to a global reduction of all voltages 
which adversely affects the accuracy of the method. This 
can be resolved by the additional introduction of weight- 
ing factors < Wj < 1 for each electrode such that the 
matrix A is replaced by A ■ ■ = AijWj. The voltages w' ob- 



tained from Eq. ( 38 1 with accordingly changed matrices 



U, V, S' and D are then to be rescaled l/wj. A reason- 
able procedure would now be to start with all Wj = 1 and 
to iteratively decrease the Wj for electrodes for which the 
voltage Uj is out of range. 



u = VS'U T ® 



VDV T u . 



(38) 



D is a N x N diagonal matrix with entries di 
which shows an opposite behavior as the Thikhonov trun- 
cation function of Eq. ( [35] ): where the truncation func- 
tion leads to vanishing voltages avoiding divergencies, di 
tends towards one (see Fig. [5] blue curve) such that the 
second term in Eq. (38) keeps the voltages close to Uq. 



By contrast for large values of Si the corresponding val- 
ues of di are vanishing as s~ 2 such that the second term 
in Eq. (35) has no effect. The N x N matrix V is needed 



http : //www. netlib . org/lapack 



B. Application 

A full implementation of the algorithm can be found 
in the supplied numerical tool box in the svdreg pack- 
age. Fig. [6] shows the obtained voltages which realize 
a harmonic trapping potential $(x) = 8(x — Xq) 2 with 
5 = 0.03y/mm 2 at different positions Xq with a given 
voltage range of —10 < Ui < 10 for the linear five- 
segmented trap of Fig. [ljb). We have also experimen- 
tally verified the accuracy of the potentials by using the 
ion as a local field probe with sub-percent agreement to 
the numerical results (Huber et al. 2010). 
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V. QUANTUM DYNAMICS - EFFICIENT NUMERICAL 
SOLUTION OF THE SCHRODINGER EQUATION 

As the trapped atoms can be cooled close to the ground 
state, their motional degrees of freedom have to be de- 
scribed quantum mechanically. In this chapter we present 
methods to solve the time-independent Schrodinger equa- 
tion and the time-dependent Schrodinger equation. The 
presented tools are used in the Sec. [VI] and |VII| about 
optimal control. For trapped ions, the relevance of treat- 
ing the quantum dynamics of the motional degrees is not 
directly obvious as the trapping potentials are extremely 
harmonic, such that (semi)classical considerations are of- 
ten fully sufficient. However, full quantum dynamical 
simulations are important for experiments outside the 
Lamb-Dicke regime ( |McDonnell et al. 2007 |Poschingerj 



et al. 2010 ), and for understandin g the sources of gate 
infidelities (Kirchmair et al. 2009 ). For trapped neutral 



atoms, the confining potentials are generally very anhar- 
monic such that quantum dynamical simulations are of 
fundamental importance. 



A. Solution of the time-independent Schrodinger equation 
- the Numerov method 

The stationary eigenstates for a given external poten- 
tial are solutions of the time-independent Schrodinger 
equation (TISE) 



2m dx 2 



<S>(x) ip(x) = Eip(x). 



(39) 



For the harmonic potentials generated with the method 
of the previous chapter these solutions are the harmonic 
oscillator eigenf unctions. But how can we obtain the 
eigenfunctions and eigenenergies for an arbitrary poten- 
tial? A typical textbook solution would be choosing a 
suitable set of basis functions and then diagonalizing the 
Hamiltonian with the help of linear algebra packages such 
as lapack to obtain the eigenenergies and the eigenfunc- 
tions as linear combination of the basis functions. 

A simple approach is exploiting the fact that phys- 
ical solutions for the wavefunction have to be normal- 
izable. This condition for the wavefunction leads to 
the constraint that the wavefunction should be zero for 
x — > ±oo. Thus it can be guessed from the potential 
shape where the nonzero parts of the wavefunction are 



located in space, and Eq. (39) can be integrated from a 



starting point outside this region with the eigenenergy as 
an initial guess parameter. For determining the correct 
energy eigenvalues, we make use of the freedom to start 
the integration from the left or from the right of this 
region of interest. Only if correct eigenenergies are cho- 
sen as initial guess, the two wavefunctions will be found 
to match (see Fig. [7|). This condition can then be ex- 
ploited by a root-finding routine to determine the proper 
eigenenergies. If the Schrodinger equation is rewri tten as 
-j-^ip(x) = g(x)tp(x), then the Numerov method (Blatt 
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FIG. 7 (Color online). Illustration of the Numerov algorithm 
for the numerical solution of the TISE: The dashed line shows 
the trapping potential for a particle. The black solid line 
shows the ground state eigenfunction. Blue and red lines show 
the result of numerical integration starting from right and left 
respectively if the energy does not correspond to an energy 
eigenvalue. 



1967) can be used for integration. Dividing the x-axis 



into discrete steps of length Ax, the wavefunction can be 
constructed using the recurrence relation 



4>n 



+ 1 



ip n (2 + ^g n Ax 2 ) - V„-i(l - ±g n -iAx 2 ) 



Ax2 n 



+0(Ax 6 



(40) 



where tp n — tp(% + nAx), g n — g{x + nAx). With this 
method, the stationary energy eigenstates are obtained. 
The source code is contained in the octtool package. 



B. Numerical evaluation of the time-dependent 
Hamiltonian 

In order to understand the behavior of quantum sys- 
tems under the influence of external control field and to 
devise strategies for their control, we perform numerical 
simulations of the time evolution. In the case of systems 
with very few degrees of freedom, the task is simply to 
solve linear first order differential equations and/or par- 
tial differential equations, depending on the representa- 
tion of the problem. Already for mesoscopic systems, 
the Hilbert space becomes so vast that one has to find 
suitable truncations to its regions which are of actual rel- 
evance. Here, we will only deal with the simple case of 
only one motional degree of freedom and possible addi- 
tional internal degrees of freedom. 

An essential prerequisite for the propagation of quan- 
tum systems in time is to evaluate the Hamiltonian; first, 
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one must find its appropriate matrix representation, and 
second, one needs to find an efficient way to describe its 
action on a given wavefunction. The first step is decisive 
for finding the eigenvalues and eigenvectors, which is of- 
ten important for a meaningful analysis of the propaga- 
tion result, and the second step will be necessary for the 
propagation itself. We assume that we are dealing with a 
particle without any internal degrees of freedom moving 
along one spatial dimension x. A further assumption is 
that the particle is to be confined to a limited portion of 
configuration space < x < L during the time interval 
of interest. We can then set up a homogeneous grid 



Xi — i Ax, i = 1, . . . , N, Ax = 



L 
N' 



(41) 



A suitable numerical representation of the wavefunction 
is given by a set of N complex numbers 



ipi = 1p(t,Xi). 



(42) 



The potential energy part of the Hamiltonian is diagonal 
in position space, and with Vi — V(xi) it is straightfor- 
wardly applied to the wavefunction: 



v^{x) -> Vii>i. 



(43) 



One might now wonder how many grid points are nec- 
essary for a faithful representation of the wavefunction. 
The answer is given by the Nyquist- Shannon sampling 
theorem. This theorem states that a band limited wave- 
form can be faithfully reconstructed from a discrete sam- 
pled representation if the sampling period is not less than 
half the period pertaining to the highest frequency oc- 
curring in the signal. Returning to our language, we can 
represent the wavefunction exactly if its energy is lim- 
ited and we use at least one grid point per antinodc. 
Of course, one still has to be careful and consider the 
possible minimum distance of antinodes for setting up a 
correct grid. Eq. (42 1 then gives an exact representation, 



and Eq. ( 43 ) becomes an equivalence. 



The kinetic energy operator, however, is not diagonal 
in position space, because the kinetic energy is given by 
the variation of the wavefunction along the spatial coor- 
dinate, i.e. its second derivative: 



T 



2m dx 2 



(44) 



One could then apply T by means of finite differences (see 
Sec. n| which turns out to be extremely inefficient as one 
wou d have to use very small grid steps (large N) in order 
to suppress errors. At the very least, we would have to be 
sure that the grid spacing is much smaller than the min- 
imum oscillation period, which is in complete contrast to 
the sampling theorem above. In order to circumvent this 
problem, we consider that T is diagonal in momentum 
representation, with the matrix elements 



T nn ' — (k n \T\k rl 



2 1.2 



h 2 k 



2m 



(45) 



Thus, we can directly apply the kinetic energy operator 
on the wavefunction in momentum space: 



i = 1, 



,M, 



(46) 



where ip{k) is the momentum representation of the wave- 
function. 

The quantity we need is the position representation of 
ip{x) with the kinetic energy operator applied to it, which 
gives 



(fV») j = (x,\t\ip) 



N 



3=1 
N M 

(xi | fc„) (k n | T\ k n ) { fc» | Xj ) (a 

j = \ n=l 



M 

M ^ 



71=1 



n 2 k 2 

2m 



N 



A? 



where J 7 , 



(k n \xj) 



(47) 



/ v M. An explicit expres- 



sion for the matrix Tij will be given below. After addition 
of the diagonal potential energy matrix one obtains the 
total Hamiltonian in the position representation, it then 
can be diagonalized by means of computer algebra pro- 
grams or efficient algorithms as the dsyevd routine of the 
computational algebra package LAPACK. 

An interesting perspective on the propagation problem 
is seen when the second last line of Eq. ( 47 1 is read from 



right to left, which gives a direct recipe for the efficient 
application of the kinetic energy operator: 

1. Transform the initial wavefunction to momentum 
space by performing the matrix multiplication 



N 



i>3 



(48) 



2. Multiply with the kinetic energy matrix elements: 

4>' n = f nn i> n . (49) 

3. Transform back to position space by performing an- 
other matrix multiplication 



(50) 



n=l 



These three steps can, of course, be merged into only 
one matrix multiplication, but the crucial point here is to 
notice that the matrix multiplications are nothing more 
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than a Discrete Fourier Transform (DFT), which can be 
performed on computers with the Fast Fourier Transform 
algorithm (FFT) flCooley and Tukey|[T965| . This has the 
tremendous advantage that instead of the N 2 scaling for 
matrix multiplication, the scaling is reduced to NlogN. 

Up to here, we have made no statement on how the 
grid in momentum space defined by the k n and M is to 
be set up. The usage of FFT algorithms strictly requires 
M = N . The grid in momentum space is then set up by 



TlAt 



N 
~2 



l<n< 



N 



Afc 



2K 



(51) 



analogously to Eq. ( 41 1 
then simply T ma: 



The maximum kinetic energy is 
1 . The remaining free parameter 



for a fixed position space grid determined by L and N is 
then the maximum wavenumbcr K. Its choice is moti- 
vated by the sampling theorem: the position space step 
Ax is to be smaller than the minimum nodal distance 
Amin/2 = ir/K, which establishes the relation 



H K AkN 



(52) 



j3 is a 'safety factor' which is to be chosen slightly smaller 
than one to guarantee the fulfillment of the sampling the- 
orem. For the sake of clarity, we note that Eq. (52) is 
equivalent to 



L 

N 



K ' 



(53) 



The optimum number of grid points iVopt for efficient 
but accurate calculations is then determined by energy 
conservation, i.e. the grid should provide a maximum 
possible kinetic energy T max equal to the maximum pos- 
sible potential energy V max . The latter is determined by 
the specific potential pertaining to the physical problem 
under consideration, whereas T max is directly given by 
the grid step. We can therefore state: 



V m 



h 2 K 2 



N, 



opt 



2m 

1 / f3irhN\ 2 
2m V L J 
L 



L 

\/2mV m 



(54) 



If more grid points are chosen, the computational ef- 
fort increases without any benefit for the accuracy. 
In turn the results become inaccurate for fewer grid 
points. If we consider a harmonic oscillator with V(x) — 

) , we obtain Kn ax = \mbj 2 L 2 and therc- 



\moj 2 (x 
fore 



N opt = 



muiL 2 

(3h : 



(55) 



which is (for f3 = 1) exactly the number of eigenstates 
sustained by the grid. 



Now, the recipe for the application of the kinetic en- 



ergy operator can be safely performed, while in Eqs. (48) 



and (50) the matrix multiplications are simply to be 
replaced by forward and backward FFTs respectively 4 . 
This Fast Fourier Grid Method was first presented by 



Feit et al. (1982) and Kosloff and Kosloff (1983). For a 



(1988) 



review on propagation schemes, see Kosloff 

With the grid having been set up, we are in good 
shape to calculate the eigenstates of a given Hamilto- 
nian by matrix diagonalization as mentioned above. The 
explicit expression for the position space matrix elements 
of Eq. (|47| is ( |Tannor[ [2007] ) 



T u = 




for I = j, 

otherwise 



(56) 



The efficiency of the Fourier method can be signifi- 
cantly increased for problems with anharmonic poten- 
tials, as is the case with the Coulomb problem in molec- 
ular systems. If one is looking at the classical trajectories 
in phase space, these will have distorted shapes, such that 
only a small fraction of phase space is actually occupied 
by the system. This leads to an inefficient usage of the 
grid, unless an inhomogeneous grid is used. The Nyquist 
theorem can be invoked locally, such that the local de 
Broglie wavelength \ dB — 2n [2m (E — V(x))] ' is to 
be used as the grid step. Further information can be 



in Refs. ( 


Fattal et al. 


Willner et al. 


2004 


)■ 



The general propagator for time-dependent Hamilto- 
nians is given by 



(57) 



with the time ordering operator T and a Hamiltonian 
consisting of a static kinetic energy part and a time- 
dependent potential energy part, which one may write 
as 



H{t) =f + V(t). 



(58) 



We discretize the problem by considering a set of interme- 
diate times i„, which are assumed to be equally spaced. 
The propagator is then given by 



U(t,t )=T~[[e-^^ At . 



(59) 



4 Two caveats for handling FFTs shall be mentioned. First, some 
FFT algorithms do not carry out the normalization explicitly, 
such that one has to multiply the resulting wavefunction with 
the correct normalization factor after the backward FFT. Second, 
one is initially often confused by the way the data returned by 
the FFT is stored: the first component typically pertains to k = 
0, with k increasing by +Afe with increasing array index. The 
negative k components are stored from end to beginning, with 
k = —Afc as the last array element, with fc changing by —Afc 
with decreasing array index. 
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It is important to state that the time-ordering operator 
now just keeps the factors ordered, with decreasing t n 
from left to right. Here, the first approximation has been 
made by replacing the control field by its piecewise con- 
stant simplification V(t n ). The short term propagators in 
the product of Eq. ( 59 ) have to be subsequently applied 



to the initial wavefunction. The remaining problem with 
the application of the short time propagators then arises 
due to the non-commutativity of T and V(t n ). A possi- 
ble way out would be the diagonalization of T + V(t n ) 
in matrix representation, which is highly inefficient due 
to the unfavorable scaling behavior of matrix diagonal- 
ization algorithms. Two main solutions for this prob- 
lem are widely used, namely the split-operator technique 
and polynomial expansion methods, which are to be ex- 
plained in the following. 



C. The split-operator method 

The basic idea of the split- operator method is to sim- 
plify the operator exponential by using the product 



,-lH(t n )At 



,-^V(t n ) At -|T At -^V(t n ) At 



(60) 



at the expense of accuracy due to violation of the non- 
comm utativity of th e kinetic and potential energy oper- 
ators (Tannor 20071. The error scales as At 3 if V(t n ) is 
taken to be the averaged potential over the time interval 
At ( |Kormann et al. 2008| ) with t n being the midpoint of 
the interval. Additional complexity arises if one is deal- 
ing with internal degrees of freedom, such that distinct 
states are coupled by the external control field. This is 
exactly the case for light-atom or light-molecule interac- 
tion processes. In these cases no diagonal representation 
of V exists in position space, meaning that it has to be 
diagonalized. 



D. The Chebyshev propagator 

A very convenient way to circumvent the problems as- 
sociated with the split-operator method is to make use 
of a polynomial expansion of the propagator 



has to be mapped on this interval by shifting and rescal- 



mg: 



H> = 2 H „- E ^-1, 



(63) 



where H is the unity matrix, and E > and E < denote 
the maximum and minimum eigenvalues of the unsealed 
Hamiltonian, respectively. The propagation scheme is 
then as follows: 

1. Given the initial wavefunction ip(t = U), set 

= i>(t = u), 

0x = -iH'<j> . (64) 



2. Calculate 



(65) 



for all n < n„ 



which is the recursion relation 



Eq. ( 62 ) applied on the wavefunction. 



3. Sum the final wavefunction according to 



= P -m( E <+ E >) At >T 



n=0 



a n 4> n . (66) 



The phase factor in front of the sum corrects for the 
energy rescaling and the expansion coefficients are given 
by Bessel functions: 



Jo 



(E > -E < )At 

2h 



for n = 0, 



K~i) n Jn( (E> -£ <)At ) fom>0. 



(67) 



exp 



It is interesting to note that the Chebyshev polynomi- 
als are not used explicitly in the scheme. Due to the 
fact that the Bessel functions J n (z) converge to zero ex- 
ponentially fast for arguments z > n, the accuracy of 
the propagation is only limited by the machine preci- 
sion as long as enough expansion coefficients are used. 
For time-dependent Hamiltonians however the accuracy 
is limited by the finite propagation time steps At, which 
should be much smaller than the time scale at which the 
control fields are changing. A detailed account on the 
accuracy of the Chebyshev method for time-dependent 



H(t') dt'j = dfc Hk(H), (61) problems is given in ( |Ndong et al] |2010| |Peskin et al. 



1994). The suitable number of expansion coefficients can 



and exploit the properties of these polynomials. As we 
will see, an expansion in terms of Chebyshev polynomials 
Hk leads to a very favorable convergence behavior and a 
simple implementation due to the recurrence relation 



Hk+i(x) = 2xH k (x) - Hk-i(x), 



(62) 



with Hq(x) = 1 and Hi(x) — x. As Chebyshev polyno- 
mials are defined on an interval x G [—1,1], the energy 



easily be found by simply plotting their magnitudes. The 
most common error source in the usage of the propaga- 
tion scheme is an incorrect normalization of the Hamilto- 
nian. One has to take into account that the eigenvalues 
of the Hamiltonian might change in the presence of a 
time-dependent control field. A good test if the scheme 
is working at all is to initialize the wavefunction in an 
eigenstate of the static Hamiltonian and then check if 
the norm is conserved and the system stays in the initial 
state upon propagation with the control field switched 
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off. Another important point is that the propagation ef- 
fort is relatively independent of the time step At. For 
larger time steps, the number of required expansion coef- 
ficients increases linearly, while the total number of steps 
decreases only. For extremely small steps, the computa- 
tional overhead of the propagation steps will lead to a 
noticeable slowdown. All presented numerical propaga- 
tors are contained in the octtool package. 



VI. OPTIMIZING WAVEPACKET MANIPULATIONS - 
OPTIMAL CONTROL THEORY (OCT) 

Now that we know how to efficiently simulate 
wavepackets in our quantum system and how to manip- 
ulate the potentials, we can begin to think about de- 
signing our potentials to produce a desired evolution of 
our wavepacket, whether for transport, or for more com- 
plex operations (such as quantum gates). In this chap- 
ter, we discuss one method for achieving this in detail, 



namely optimal control theory 


Khaneja et al. 2005 


Kosloff et al. 1989 Krotov 2008 


1996 Peirce et al. 


1988 Sklarz and Tannor 2002 Som 


oi et al. 1993 Tan- 


nor et al. 1992 Zhu and Rabitz 


1998). These methods 



belong to a class of control known as open-loop, which 
means that we specify everything about our experiment 
beforehand in our simulation, and then apply the results 
of optimal control directly in our experiment. This has 
the advantage that we should not need to acquire con- 
stant feedback from the experiment as it is running (an 
often destructive operation in quantum mechanics). We 
will focus on one particular method prevalent in the lit- 
erature, known as the Krotov algorithm ( |Krotov| |1996 
Somloi et al] 119931 iTannor et oZ.l 119921). 



A. Krotov algorithm 

Optimal Control Theory (OCT) came about as an ex- 
tension of the classical calculus of variations subject to 
a differential equation constraint. Techniques for solv- 
ing such problems were already known in the engineer- 
ing community for some years, but using OCT to op- 
timize quantum mechanical systems only began in the 
late 1980s with the work of Rabitz and coworkers ( |Peirce 
et al. 1988 Shi et al. |1988 ), where they applied these 
techniques to numerically obtain optimal pulses for driv- 
ing a quantum system towards a given goal. At this 
time, the numerical approach for solving the resulting 
set of coupled differential equations relied on a simple 
gradient method with line search. In the years that fol- 
lowed, the field was greatly expanded by the addition of 
more sophisticated techniques that promised improved 
optimization performance. One of the most prominent 
amongst these is the Krotov method ( Krotov 1996 Som- 



loi et al. 1993 Tannor et al. 1992) developed by Tan- 



nor and coworkers for problems in quantum chemistry 
around the beginning of the 1990s, based on Krotov's 
initial work. This method enjoyed much success, being 



further modified by Rabitz ( Zhu and Rabitz 1998 ) in the 
late 90s. 

Until this point, OCT had been applied mainly to 
problems in quantum chemistry, which typically involved 
driving a quantum state to a particular goal state (known 
as state-to-state control), or maximizing the expectation 
value of an operator. The advent of quantum information 
theory at the beginning of the new millennium presented 
new challenges for control theory, in particular the need 
to perform quantum gates, which are not state-to-state 
transfers, but rather full unitary operations that map 
whole sets of quantum states into the desired final states 
dCalarco et al. 1 120041 iHsieh and Rabitzl 120081 iPalao and 



Kosloff 2002). Palao and Kosloffl (|2002Textended the 



Krotov algorithm to deal with such unitary evolutions, 
showing how the method can be generalized to deal with 
arbitrary numbers of states. This will become useful for 
us later when we want to optimize the Cirac-Zoller gate. 
Other methods for optimal control besides Krotov have 
also been extensively studied in the literature, mos t no- 
tably perhaps being GRAPE (Khaneja et al. 2005), but 



these methods will not be discussed here. 



a. Constructing the optimization objective We will now 
proceed to outline the basics of optimal control theory 
as it is used for the optimization later in this paper 5 . 
We always begin by defining the objective which is a 
mathematical description of the final outcome we want 
to achieve. For simplicity, we shall take as our example a 
state-to-state transfer of a quantum state \ip(t)} over the 
interval t S [0, T]. We begin with the initial state \i/j(0)), 
and the evolution of this state takes place in accordance 
with the Schrodinger equation 



(68) 



where H is the Hamiltonian. (Note that we will often 
omit explicit variable dependence for brevity.) Now as- 
sume that the Hamiltonian can be written as 



H(t) = H {t) +J2 £ i{t)Mt), 



(69) 



where Hq is the uncontrollable part of the Hamiltonian 
(meaning physically the part we cannot alter in the lab) , 
and the remaining Hi are the controllable parts, in that 
we may affect their influence through the (real) func- 
tions £i(t), which we refer to interchangeably as 'controls' 
or 'pulses' (the latter originating from the early days of 
chemical control where interaction with the system was 
performed with laser pulses). Let's take as our goal that 
we should steer the initial state into a particular final 



For a tutorial on quantum optimal control theory which covers 
the top ics presented here in more detail, see Wcrschnik and Gross 
(2007i. 
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state at our final time T, which we call the goal state 
\<j>). A measure of how well we have achieved the final 
state is given by the fidelity 



(70) 



which can be seen simply as the square of the inner prod- 
uct between the goal state and the final evolved state. 
Note that J\ [ip] is a functional of ip. 

The only other constraint to consider in our problem 
is the dynamical one provided by the time-dependent 
Schrodinger equation (TDSE) in Eq. (68). We require 



that the quantum state must satisfy this equation at all 
times, otherwise the result is clearly non-physical. If a 



quantum state \t/>(t)) satisfies Eq. (68 1, then we must 
have 

(d t + i£) m)) =0,Vte T, where d t = ^. (71) 

We can introduce a Lagrange multiplier to cast our con- 
strained optimization into an unconstrained one. Here, 
we introduce the state \x(t)) to play the role of our La- 
grange multiplier, and hence we write our constraint for 
the TDSE as 



J2hA,x}= f (<x(*)l($ + £#M*)) + c.c.) 



at 



= 2Re / ( X (t)\(dt + iH)\iP(t)) dt, (72) 
Jo 

where we have imposed that both \ip(t)} and (tp(t)\ must 
satisfy the TDSE. 



b. Minimizing the objective Now that we have defined 
our goal and the constraints, we can write our objective 
J[ei,ip,x] as 



(73) 



The goal for the optimization is to find the minimum of 
this functional with respect to the parameters ip{t), x(t) 
and the controls £i(t). In order to find the minimum, 
we consider the stationary points of the functional J by 
setting the total variation 5 J = 0. The total variation is 
simply given by the sum of the variations SJ^ (variation 
with respect to (variation with respect to %), and 

8J 6i (variation with respect to £;), which we set individ- 
ually to zero. For our purposes, we define the variation 
of a functional 



S^F[ij;} = F[ip + Sip] - F[ip]. 



(74) 



This can be thought of as the change brought about in 
F by perturbing the function ip by a small amount Sip. 
Considering SJ^p, we have 



Using our definition from Eq. (74 1 for 5J\^ results in 
SJirt = Ji[ip + 5ip] - Ji[ip] 

= -|(0|^r) + ^(r))| 2 + |(^(T))| 2 

= -(ip(T) + 8^{T)\cP){^{T) + 5^{T)) 
+ I«#(T)>| 2 

= -<ty(T)|M^(T)> - MTMM^CT)) 
-\(cf>\SiP(T))\ 2 

= -2Re{(^(T)|0)(<^(T))} - \(m(T))\ 2 - 

(75) 

The last term is 0(Sip(T) 2 ), and since 8ip(T) is small we 
set these terms to zero. Hence, we have finally 



<yJi^ = -2Re{<V(T)|^)^|^(T))}- 
Repeating this treatment for SJ2,4>, we have 



(76) 



SJ^ = 2Re / ( X {t)\(d t + lH)\ip{t) + S^(t))dt 
Jo 

-2Re [ ( X (t)\(dt + iH)m))dt 
Jo 

= 2Re / ( X (t)\(dt + jiH)\Sm)dt 
Jo 

= 2Rc{(x(T)|^(T))-( X (0)|^(0)) 

(( x (t)\iH - (d tX (t)\) \Sij(t))dt}. 

(77) 

Noting that the initial state is fixed, we must have 
Sip(0) = 0. Thus setting 5J^ — 0, we obtain the two 
equations 

(^(T)|0)(^(T)) + (x(T)\Sip(T)} = 0, (78) 
(( X (t)\iH-(d t x(t)\)\SiP(t)}=0. (79) 

Since these must be valid for an arbitrary choice of \Sip), 
we obtain from Eq. ( 78 ) the boundary condition 



\x(t)} = mmn, 



and from Eq.(79) the equation of motion 

ihd t \ x (t))=H\x(t)). 



(80) 



(81) 



5 Jih = 5Ji^ + 5 J: 



2,ip- 



We now continue the derivation by finding the vari- 
ation S J X , which results in the condition already given 
in Eq. (68), namely that \ip) must obey the Schrodinger 
equation. The variation SJ £i results in the condition 



llm \m\\^m)) =0, (82) 



where can be used to suppress updates at t = and 
t = T. It can be clearly seen, that Eq.(82) cannot be 
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solved directly for the controls e< since we have a system 
of split boundary conditions: \ip(t)} is only specified at 
t = 0, and similarly |x(*)) only at t = T. Hence we 
require an iterative scheme which will solve the equations 
seli-consistently. 



c. Deriving an iterative scheme The goal of any iterative 
method will be to reduce the objective J at each iteration 
while satisfying the constraints. Written mathematically, 
we simply require that J k+1 — J k < 0, where J k is the 
value of the functional J evaluated at the fcth iteration of 
the algorithm. We will also attach this notation to other 
objects to denote which iteration of the algorithm we are 
referring to. Looking at Eq. (82) and taking into account 



our constraints, the optimization algorithm presents itself 
as follows: 

1. Make an initial guess for the control fields £i(t). 

2. At the fcth iteration, propagate the initial state 
IV'(O)) until \ip k (T)) and store it at each time step 
t. 



3. Calculate \x k {T)) from Eq. (80b. Our initial guess 



from step (1) should have been good enough that 
the final overlap of the wavefunction with the goal 
state is not zero (otherwise Eq. ( 80 ) would give us 
\x k (T))=0). 

4. Propagat e \x k {T)) backward in time in accordance 
with Eq. (81 ). At each time step, calculate the new 



control at time t 



Jfe+i 



(t)=s k (t) +1 -lm 

A 7 ; 



(X k (t)\l^ k (t))\, (83) 



5. Start again with | -0(0)} , and calculate the new con- 
trol at time t 



-k+l 



(t)=e k {t)- 



A, 



Im 



(x'Wligl^'W) 



(84) 



Use these new controls to propagate |-0 fc+1 (O)) to 
obtain \*P k+1 (T)). 

6. Repeat steps (3) to (5) until the desired conver- 
gence has been achieved. 

This new method looks very similar to the gradient 
method, except that now we see that we must not use 
the 'old' \ip k (t)) in the update, but the 'new' \ifj k+1 (t)). 
This is achieved by immediately propagating the cur- 
rent \ip k+1 (t)) with the newly updated pulse, and not 
the old one. To make this explicit, take at t = 0, 
l^+^O)) = 1^(0)). We use this to calculate the first 
update to the controls s k+1 (0) from Eq. (|84|). We use 



these controls to find \ip k+1 (At)) , where At is one time- 
step of our simulation. We then obtain the next update 



e k+1 (At) by again using Eq. (84 1 where we use the old 



\x k {At)} that we had saved from the previous step. In 
other words, |x(i)) is always propagated with the old 
controls, and \ip(t)} is always propagated with the new 
controls. 

For a full treatment of th is method in the literature, see 
Sklarz and Tannor (2002). For our purposes, we simply 



note that the method is proven to be convergent, mean- 
ing J k+1 — J k < 0, and that it has a fast convergence 
when compared to many other optimization algorithms, 
notably the gradient method (Somloi et al. 1993). 



B. Application 



which one can identify as a gradient-type algo- 
rithm, using Eq. (82) as the gradient. The param- 
eter 7 is determined by a line search to ensure our 
condition that J k+1 - J k < 0. 

5. Repeat steps (2) to (4) until the desired conver- 
gence has been achieved. 

This gradient-type method, while guaranteeing conver- 
gence, is rather slow. A much faster method is what is 
known as the Krotov method in the literature. Here, the 
modified procedure is as follows: 

1. Make an initial guess for the control fields £i(t). 

2. Propagate the initial state \ip{0)) until \ip(T)). 

3. At the fcth iteration, calculate |x fc (T)) from 
Eq. ( |80"| (again taking care that the final state over- 
lap with the goal state is non-zero). 

4. Propagate \x k (T)) backward in time in accordance 
with Eq. (81) to obtain |x(0)), storing it at each 
time step. 



We now show how the wavefunction can be controlled 
via tailored time-dependent external potentials. In ion- 
trap experiments quantum control of the wavefunction 
via individual electrodes is difficult, as these are several 
orders of magnitude larger than the size of the wave- 
function. However, matching length scales occur for cold 



atoms trapped in optical lattice potentials (Calarco et al. 



2004) or magnetic micro trap. We nonetheless focus on 



the ion trap system and present this as a generic exam- 
ple for quantum control where we would like to transport 
the ion from one place to another without exciting higher 
motional states. The transport is performed by applying 
time-dependent voltages Ui(t) to the electrodes generat- 
ing a total potential 

n 

$(x, utit), u 2 (t), u n {t)) = $i{x)Ui(t). (85) 

1=1 

The Hamiltonian of this system is 

ft 2 d 2 



H {x,ui{t), ...,u n (t)) 



2m dx 2 



$(x,iti(f),...,u n (t)). 

(86) 
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qubit and are manipulated by laser or microwave radi- 
ation. They are subject to quantum dynamics where a 
certain degree of sophistication can arise due to the pres- 
ence of interference effects. Joint manipulation of the 
internal and external degrees of freedom is therefore a 
promising playground for the application of OCT theory 
along with the quantum dynamical simulation methods 
presented in chapter |Vj In the following we present the 
Cirac-Zoller controlled-not gate as a case study. We first 
explain the gate mechanism, then we use the OCT al- 
gorithm to derive a laser pulse sequence for the proper 
realization of the quantum gate. 



A. The Cirac-Zoller controlled-not gate 



FIG. 8 (Color online). Fidelity increase after iterative opti- 
mization steps, (a) Initial guess of the time-dependent volt- 
ages applied to the electrodes, (b) Resulting wavefunction at 
the final time T. The fidelity is only 0.3. (c) Voltage configu- 
rations obtained after 100 iterative calls of the Krotov optimal 
control method, (d) The final wavefunction obtained at time 
T with optimized voltages agrees well with the target ground 
state wavefunction. The fidelity is increased to 0.997. 



As a target wavefunction \4>), we choose a harmonic os- 
cillator ground state wavefunction centered at the target 
position. We thus have to maximize the wavefunction 
overlap \(<fi\ip(T))\ 2 , where \ip(T)) is the wavefunction 
at the final time T obtained by application of the time- 
dependent voltages. This is exactly the fidelity functional 
of Eq. ( 70 1 and the initial condition for the Langrange 



multiplier of Eq. ( 80 ) . The update Eq. ( 84 ) for the con- 
trol parameters ufJT) from iteration step k to k + 1 has 
the form 



,/c+l 



(*) = u k S) + ^{x k {t)\- h Ux)\^ k+1 {t)). (87) 



Starting with a sinusoidal-shaped initial guess for 
the time-dependent voltage configuration u®(t) (see 
Fig. |8][a)) the wavefunction at the final time T is ex- 
cited to the first excited state, leading to a wavefunc- 
tion overlap of 0.3 as seen in Fig. |8jb). After 100 steps 
of optimization the wavefunction overlap has been itera- 
tively increased to 0.997 (Fig.[8| and the motional ground 
state has been preserved (see Fig. [8](d)) . The optimized 
time-dependent voltages can be seen in Fig. [8jc) . The 
full source code of the optimal control algorithm to per- 
form the presented optimization is contained in the oct- 
tool package. 



VII. IMPROVING QUANTUM GATES - OCT OF A 
UNITARY TRANSFORMATION 

Up to now we have only considered the motional de- 
grees of freedom of trapped ions. The internal electronic 
degrees of freedom serve as information storage for the 



Atomic qubits are suitable carriers of quantum infor- 
mation because one can exploit the internal electronic 
degrees of freedom for realizing a near-perfect two level 
system ( |Cohen-Tannoudji et al. 2004), where coherent 
transitions between the states can be driven by laser radi- 
ation. Suitable long lived states are found to be sublevels 
of the electronic ground state connected via stimulated 



Raman transitions ( Poschinger et al. 2009 ) or metastable 



electronic states excited on narrow dipole forbidden tran 



sitions ( Schmidt-Kaler et al. 2003a I. For neutral atoms 



the additional possibility of employing Rydberg states ex- 
ists ( Gaetan et al. 2009 ). In the following we will explain 
the Cirac and Zoller ( 1995 ) controlled-no t (cnot) gate as 

(|2003b|). To understand 



realized bylSchmidt 



<aler et al. 



each single step of the gate operation we first h ave to be- 
comc acquainted with the light- ion interaction (Lcibfried 



et al. 2003). In the following \l) and |f) denote the qubit 
states with energies and fku , respectively. The full 
Hamiltonian of the system is H = Hq + H a + Hl , where 



Hq is given by Eq. ( 86 ) with a constant harmonic trap 
potential $(x) — muj^x 2 /2, H a = |f) (t| f^o describes 
the energy of the internal electronic excitation. de- 



scribes the interaction betwee n light and atom (James 
" 2002h 6 : 



1998 



Sasura and Buzck 



(It) (11 + 11) (t|) e 



Hr = 



where f2 is the Rabi frequency of the transition between 
the qubit states. In this form, the Hamiltonian contains 



6 For a dipole transition can be written as H'g = — er ■ E, and 
similarly for a quadrupole transition = — | . r^Tj , 
with E = Erjcos(k • x — w^t — </>), where r denotes the relative 
position of the valence electron with respect to the nucleus and 
x is the position of the ion. The frequency of the laser is given 
by u>l, with the optical phase </>■ We obtain the matrix elements 
of Hl by means of the identity operator 11 = |4) (4,| + |t) such 
that Hl = llff^U. All diagonal elements vanish and the Hamil- 
tonian can be written as in Eq. where we have expressed 
the cosine by exponentials. 
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terms oscillating at the laser frequency. For optical fre- 
quencies we would need rather small time steps for an 
accurate numerical simulation. To avoid this, we change 
to the interaction picture 



H, 



e lH « t/h \iJj), and 
e iH a t/h He -iH a t/n 



(89) 
(90) 



where is the state of the motional and internal degrees 
of freedom in the Schrodinger picture. Hj can be ex- 
panded by using the Baker-Campbell-Hausdorff formula 7 



Ht — Hn 



o iu t 



|t> (|| +h.c.) f e ^ kx - ULt - 



h.c. 



(. 91 ) 

If we additionally make the rotating wave approximation, 
i.e. neglect fast oscillating terms at the frequency lul+ljq, 
we obtain 



Hi = H 



Ml 



It) (II 



a i(kx—5t- 



h.c. 



(92) 



with the laser detuning S = wl — ujq. Now we can do 
the numerics with much larger time steps. The term 
proportional to ||) (|| e lkx is responsible for absorption 
processes: it changes the ||) state into the ||) state and 
displaces the motional state in momentum space by the 
photon recoil hk. The Hermitian conjugate term propor- 



tional to ||) (t| ' 



-ikx 



is responsible for stimulated emis- 



sion processes from the ||) state back to the ||) ground 
state, where a photon is emitted back into the laser field 
displacing the motional state in momentum space by 
-Hk. 

If the laser frequency wl is tuned to the atomic res- 
onance, such that 5 = 0, then the interaction describes 
simple Rabi oscillations between the qubit states with 
frequency tt. No net momentum is transferred as ab- 
sorption processes and stimulated emission contribute 
equally. Thus, we can use this interaction for direct con- 
trol of the internal states without changing the motional 
state. If the laser is irradiated on the ion during the 
time t such that fit = tt/2 (referred to as a ir/2 pulse), 
then superposition states are created 8 : ||) — > ||) + |t) an d 
|t) -> -||) + It)- These states evolve as ||) + e - l " ot |t)- If 
we now apply a second tt/2 pulse in phase with the oscil- 
lating superposition, we obtain the states ||) + It) ~ Ht) 
and — 1|) + ||) — > — 1|). If the optical phase is shifted by 
7r such that the laser field and the superposition are oscil- 
lating out of phase, we reverse the action of the first tt/2 
pulse: ID + lt) -> ||) and -||) + |t) -> It)- Hence, we ob- 
tain orthogonal results depending on the phase relation 
between the laser and the superposition state, which is 
the basic principle of Ramsey spectroscopy. If the laser 



frequency is kept perfectly resonant and phase stability 
is maintained, then one can detect externally induced 
phase flips of the superposition state during the waiting 
time T by mapping the phase to the two states ||) and 
||). Starting from the ground state, application of reso- 
nant tt/2 pulses continuously cycles through the series of 
states 



ID 

tt/2. 



tt/2 



II) + It) 



tt/2 



-ID + lt)^ 



It) 
ID- 



(93) 



We can see that a concatenation of four it/ 2 pulses (a 2tt 
pulse) takes the system back to the ground state, however 
with a phase flip of 7r, which is due to the fundamental 
4-7T rotational symmetry of spin-1/2 systems 9 . 

When the laser frequency is tuned below the atomic 
resonance by the vibrational trap frequency w t r, i-e. the 
laser is red detuned by 5 = — u; tr , we drive red-sideband 
transitions between the states ||, n) and ||,n— 1) with 
reduced Rabi frequency flr/y'n. Rabi oscillations are ob- 



tained in a similar manner as in Eq. ( 93 1 by replacing 
||) -> | |,n) and |f) ||,n-l). n = 0,1,2,... de- 
notes the harmonic oscillator eigenstates, and r\ = xq ■ k 
is the Lamb-Dicke parameter, where xo is the size of the 
ground state wavefunction. This parameter sets the cou- 
pling strength between the laser radiation and the atomic 
motion. Atomic excitation on the red sideband is accom- 
panied by the lowering of the harmonic oscillator energy 
by one vibrational quantum. 

For a blue laser-detuning with S = w tr , the blue side- 
band interaction is realized. On the blue sideband tran- 
sitions between the states ||, n) and ||, n + 1) are driven 
with Rabi frequency flr/^/n + 1. When applying a ir- 
pulse to the ||, 0) state we excite one motional quantum 
and obtain the state ||, 1). On the other hand, if applied 
to the |t, 0) state, no motional quantum can be excited 
thus this state does not couple to the blue sideband. A 
7r pulse on the blue sideband transition can therefore be 
used to map quantum information back and forth be- 
tween the internal state of a specific ion and the motional 
state of an ion chain if a collective vibrational mode is 
driven. This operation is referred in the following as a 
SWAP operation. 

Now we have all the tools at hand to put the quantum 
gate together. In the following one ion will be referred to 
as control ion whose internal state is denoted by | • ) c , and 
a second target ion with internal state | • ) t . We want to 
flip the state of the target ion conditionally to the state 



7 e A Be~ A = Y^[A,B] m l/m\ with [A, B] m = [X,[X,Y] m -l], 
the commutator [A, B] = AB — BA and [A, B] Q = B. 

8 We have subsequently omitted all normalization factors 1 / \/2 as 
they do not change the physical interpretation. 



A global phase does not have any physical significance and can al- 
ways be absorbed into the definition of the states \ip)' —> e'"f'\ip). 
The absolute laser phase does not matter before the first laser 
pulse starts, of importance is the relative phase of the superpo- 
sition (imprinted by the first laser pulse) and a subsequent laser 
pulse. 
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of the control ion, realizing the CNOT truth table: 



ll)c|t)t 
lt}c|t}t 
It)c|t>t 



ll)c|t)t, 

lt)c|t)t, 
It)c|t>f 



(94) 



This gate is performed by the following steps: First the 
state of the control ion is mapped on a collective vibra- 
tional mode of the two ions by means of a SWAP oper- 
ation. The remaining task is to perform a CNOT gate 
between the vibrational mode and the target ion and fi- 
nally perform the SWAP -1 operation to restore the state 
of the control ion (see Fig. 9][a)). The CNOT gate between 
the motional mode and the internal state of the target 
ion corresponds to the truth table 



lt,o) t 
M)t 
lt,i)t 



|t,o> t , 
lt,o) t , 
M)t, 
It, iK, 



(95) 



where the motional state now acts as the control. The 
key element of this operation is a controlled phase gate 
between these two degrees of freedom, which corresponds 
to the truth table 



M) t 
lt,o) t 
ll,i)t 
lt,i)t 



-M) t , 

|t,0)t, 

-M)t, 
-It,i)t- 



(96) 



As explained above, mapping of a superposition phase 
to the states is accomplished by means of resonant 7r/2 
pulses. If the controlled phase is sandwiched between 
two such pulses, then the internal state of the target ion 
is given by the conditional phase from the phase gate 
operation. This realizes the CNOT gate of Eq. ( 95 ) . 



The phase gate itself is performed by exploiting on the 
one hand the fact that the blue sideband does not couple 
to the |t, 0) state and on the other hand that a 2tt pulse 
flips the phase of any given state by tt as can be seen 
from Eq. (J93b. Therefore a 2n pulse on the blue side 
band changes the phase for the states ||, 0), |t, 1) and 
|t, 1), whereas |f, 0) is left unchanged. This would result 
in the conditional phase gate of Eq. ( 96 ) such that the 



whole CNOT gate sequence is complete. 

Additional complications arise due to the fact that the 
blue sideband Rabi frequency on the |t, 1) — > It, 2) state 
is larger than the one on the |t, 0) — s- |t, 1) transition 
by a factor of \[2. The problem was resolved by a the- 
oretical proposal by |Childs and Chuang (|2000|) and re- 
alized experimentally by |Schmidt-Kaler et al.\ ( |2003b|c[ ) 
through the application of a composite pulse sequence of 
blue sideband pulses with different durations and phases 
as seen in Fig. [9]jb). In the next section we will demon- 
strate how such sequences can be automatically obtained 
by application of quantum optimal control techniques. 



(a) quantum circuit: 
ionl ||),|T) - 



SWAP 




SWAP" 1 






— i 


i — 



ion 



2 |i),|T> 



control bit 



l°> 



target bit 



ion 1 



ion 2 



blue ■* 


pulse duration/intensity 

, optical phase 


blue 


— 




7T 

7T 



t/4 



blue 

-t/4 



blue 

tt/4 



blue 

-JT I A 



FIG. 9 (a) Quantum circuit for a CNOT gate between two ions 
realized by a SWAP operation on the control ion, a CNOT gate 
between the motional state and the target ion and a SWAP - 
operation, (b) Composite pu lse sequence to realize the tota l 
CNOT gate between two ions ( Schmidt-Kaler et al. 2003b|c I. 



B. Krotov algorithm on unitary transformations 

Now we show how the optimal control algorithm from 
Sec. |VI| finds a control sequence which realizes the con- 
trolled phase gate. The input for the optimal control al- 
gorithm is the system dynamics governed by the TDSE 
and the Hamiltonian Hj from Eq. (92). The subject to 
be controlled is the unitary transform of Eq. ( 96 ) . The 



state of the system is represented by two distinct wave- 
functions in position space for the states |f) and |t)- We 
now perform the optimal control algorithm for a uni- 



tary transformation along the lines of Palao and Kosloff 
(2002). Instead of one initial and one target state, we 



now have four initial states \ip s (0)) and four target states 
\<p s ) (s — 1...4) corresponding to the states in Eq. (96). 



Additionally, we have to change the fidelity objective of 
Eq. ( 70 ) to a phase sensitive one with 



X>.ltf.(r)) 

vs=l / 



1 

2' 



(97) 



such that we have Ji[ip] = —F[i/j]. The constraint is now 
that the TDSE is to be fulfilled for all four states, thus 



we introduce four Lagrange multipliers \Xs{t))- Eq. (72) 
is now changed into 



J 2 [ £l ,^, X ] = 2^Re / (xsitm + ^H^sit^dt 

s=l J 

(98) 



From these equations wc derive the initial condition for 
the Lagrange multipliers (% S (T)| = (^> s | and the update 
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equation 



A,; 




{x k M\{^ k s + \t)) 



(99) 



C. Application 



We now compare the four pulses on the blue sideband 
(see Fig. |)Jb)) to the pulse found by the optimal con- 
trol algorithm. For our gate optimization problem we 
need only one control parameter which is the phase <f> 
of the laser on the blue sideband £\(t) = <fi(t) with the 
initial guess <j>(t) = 0. 



Fig 10 



shows the increase of 
the fidelity from 0.43 to 0.975 after 50 iterations. The 



composite pulse sequence used in ( |Schmidt-Kaler et al. 



2003b|. 



achieves a fidelity of 0.994. In both cases, the 
deviation from unity is caused by off-resonant excitation 
on the carrier transition. It is quite remarkable that the 
Krotov algorithm finds a time-dependent phase 4>(t) with 
similar amplitude and period as the composite pulse se- 
quence (see Fig. [Toj^b)). We have illustrated the OCT 
method on an example of a quantum control problem 
where a solution is already known. For more complex 
control problems however, this might not be the case, 
such that OCT allows for tackling control problems where 
the vastness of Hilbert space and the complexity of quan- 
tum interference are hard to handle manually. All sources 
can be found in the octtool package, where additionally 
a simulation of the gate operation on the wavcfunction is 
visualized. Further application of optimal control tech- 
niques to quantum gates with ions can be found in th e 
literature (IGarci'a-Ripoll et aLl 120031 IZhu et aU [2006a). 



VIII. CONCLUSION 

Applicability to other qubit types: Currently, trapped 
ion quantum systems are leading experimental efforts in 
quantum information theory. But with the growing ma- 
turity of quantum information experiments with trapped 
neutral atoms or solid state systems, we stress that the 
methods presented here will be applicable also to these 
systems with minor modifications. In this section we 
briefly elucidate how far each of the methods presented is 
applicable to each type of qubit and indicate, where ap- 
propriate, the connections between them with relevant 
citations. The methods from Sec. IIHI for the numerical 
calculation of particle trajectories are directly applicable 
to neutral atoms, where it might also be interesting to an- 
alyze their motional b ehavior in trap structures like mag - 
netic microchip traps (Fortagh and Zimmermannl 20071, 



which have grown greatly in complexity, essentially real- 
izing labs on atom chips. The ion trap community has 
mimicked the success story of atom chips by the devel- 
opment of microstructured ion traps. For neutral atoms, 
the full quantum dynamical simulations are of an even 




20 30 

iterations 

FIG. 10 (Color online) (a) Increase in fidelity from 0.43 to 
0.975 after 50 iterations (black). The composite pulse se- 
quence realizes a fidelity of 0.994 (red), (b) Time-dependent 
phase applied on the blue sideband. The dotted curve shows 
the inital guess <j>(t) = and the black curve shows the result 
obtained after 50 iterations of the Krotov algorithm. Note 
the similarity in the obtained phases concerning amplitude 
and period when compared to the compo site pulse sequence 
as used in ( Schmidt-Kaler et al. 2003b|c I (red curve). 



higher importance than for ion trap systems, which is 
due to the fact that neutral atoms are generally more 
weakly confined and are therefore much more sensitive 
to anharmonicities of the trap potentials. Quantum dy- 
namical simulations have been used in conjunction with 
the optimal control method from Sec. |VI| to investigate 
the perspectives for transport and splitting operations in 



magnetic microtraps and optical lattices ( Calarco et al. 



2004 


Chiara et al. 2008 Hohenester et al. 2007 Riedel 


et al. 


2010 


Treutlein et al. 


2006 


). The optimal con- 



trol method might turn out to be of essential importance 
for the realization of robust high-fidelity quantum gates 
in artificial atom systems like Josephson junction-based 
qubits or impurity-based qubits in solid state host ma- 
terials (Kane 19981 INeumann et al. 2008), where the 



level scheme, the coupling to external control fields and 
the decoherence mechanisms are generally more compli- 



cated thaiifortrapped ion qubits (Montangero et al. 
2007 Sporl et al. 2007). In general, the ability to cal- 
culate the quantum dynamics for any kind of engineered 
quantum system is of fundamental importance, and the 
methods presented in Sec. [V] can be directly adapted to 
any given Hamiltonian. The methods for the precise and 
efficient calculation of electrostatic fields from Sec. III. Al 
are also of interest for the optimization of electrode ge- 
ometries for quantum dot qubits based on a two dimen- 
sional electron gas (Koppens et al. 2006). Furthermore, 



the Laplace equation is of a similar mathematical struc- 
ture as the Helmholtz equation for ac electromagnetic 
fields, such that these fields may be calculated in minia- 
turized microwave traps for neutral atoms or microstruc- 
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tured Josephson devices based on adapted fast multipole 



methods (Gumerov and Duraiswami, 2004) 



Summary: We have presented the whole process of 
ion trap quantum computing starting from trap design, 
trapping, and transport of ions, and ending with laser- 
ion interactions and the simulation and optimization of 
the Cirac-Zoller CNOT gate starting from basic principles. 
The explanation of the physics is complemented by a 
detailed description of the numerical methods needed to 
perform precise simulations, which are carefully selected 
such that both precision and efficiency are maintained. 
Additionally, the numerical methods are presented in a 
general manner, such that they may be applied to solve 
physical problems outside of the particular focus of this 
paper. The source code of all methods together with all 
needed libraries packed into a single installation file can 
be downloaded from http : / /kilian-singer . de/ ent for 
both Linux and Windows operating systems. For fast 
testing of the routines we have bundled them with the 
C++ scripting library ROOT 10 . By making these libraries 
accessible to the public we want to inspire students to 
experiment, and provide researchers with a foundation 
to perform more sophisticated simulations. 
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